docs/rxd-tutorials/Restricting reactions to certain sections.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Restricting reactions to a part of a region"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A version of this notebook may be run online via Google Colab at https://tinyurl.com/neuron-restricted-reaction\n",
    " (make a copy or open in playground mode)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "Reactions may be restricted to only occur in a given part of an `rxd.Region` by multiplying the rate by an indicator variable, that is an `rxd.Parameter` that takes the value `1` in the subregion where the reaction is to occur and `0` elsewhere.\n",
    "\n",
    "Let's consider an example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup NEURON library and imports"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's import our usual NEURON libraries and definitions. Remember you can use either `um` or `µm` for micron.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from neuron import h, rxd\n",
    "from neuron.units import mV, ms, um, mM\n",
    "\n",
    "## needed for standard run system\n",
    "h.load_file(\"stdrun.hoc\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now import `plotly`, a graphics library. (You could easily modify this code to use other graphics libraries like `matplotlib`, `plotnine`, or `bokeh`.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import plotly.graph_objects as go"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Morphology and discretization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "left = h.Section(name=\"left\")\n",
    "right = h.Section(name=\"right\")\n",
    "\n",
    "left.nseg = right.nseg = 101\n",
    "left.L = right.L = 101\n",
    "\n",
    "right.connect(left)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The rxd.Region"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our region will be present on both the `left` and the `right` sections. Since they define our entire \"cell\", we will pick one and put the `rxd.Region` on the `.wholetree()` that is connected to it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cytosol = rxd.Region(left.wholetree())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The rxd.Species"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's begin by defining our initialization rule. The concentration will initially be 0 everywhere except for nodes between 90 and 110 microns from the left edge of `left` (that is, the last 10 microns of `left` and the first 10 microns of `right`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def initial_protein(node):\n",
    "    if 90 * um < h.distance(node, left(0)) < 110 * um:\n",
    "        return 1 * mM\n",
    "    else:\n",
    "        return 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the species definition itself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "protein = rxd.Species(cytosol, d=1, initial=initial_protein)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The rxd.Reaction and its localization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need an `rxd.Parameter` indicator variable to indicate the sections where the reaction is to take place. In parallel to how we defined the initial concentration of the species above, it is often simplest to define the values using a function of the node (here 1 in the `right` section; 0 otherwise):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def active_region_value(node):\n",
    "    if node in right:\n",
    "        return 1\n",
    "    else:\n",
    "        return 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the indicator variable `rxd.Parameter` itself:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "in_region = rxd.Parameter(cytosol, value=active_region_value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the reaction (here an `rxd.Rate` but could be an `rxd.Reaction` or `rxd.MultiCompartmentReaction` as well), note that we multiply by the value of `in_region`, so when that value is `1` the reaction happens and when it is `0`, the reaction does not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "production_rate = 0.002 * mM / ms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reaction = rxd.Rate(protein, production_rate * in_region)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run and visualize the simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Always initialize your simulations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.finitialize(-65 * mV)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll plot the protein concentrations against their distance from the left edge of the left section (i.e. from `left(0)`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot(fig):\n",
    "    y = protein.nodes.concentration\n",
    "    x = [h.distance(node, left(0)) for node in protein.nodes]\n",
    "    fig.add_trace(go.Scatter(x=x, y=y, name=f\"t = {h.t:g} ms\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the initial conditions, then advance for 25 ms, plot, advance again, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = go.Figure()\n",
    "plot(fig)\n",
    "\n",
    "for advance_count in range(1, 5):\n",
    "    h.continuerun(advance_count * 25 * ms)\n",
    "    plot(fig)\n",
    "\n",
    "fig.update_layout(xaxis_title=\"Position (µm)\", yaxis_title=\"Concentration (mM)\")\n",
    "\n",
    "fig.show(renderer=\"notebook_connected\")"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}